

function [obj,resL,gamma_hat]=MM_Exp(q, p, y, x, toler, gamma)

  % Returns the QR estimator for a quantile regression model under an
  % exponential specification
  % The QR estimator is computed using a MM (Majorize-Minimize) algorithm
  % from Hunter Lange (2000) JCGS
 
  %******************* Input Variables **************************************

  % q : quantile level
  % p : dimension of the covariate (including intercept)
  % y : dependent variable
  % x : covariate (must include intercept)
  % toler : tolerance parameter for the precision of the estimation
  % gamma : initial set of value for the intercept and slope coefficients
 
  %******************* Output Variables **************************************
 
  % obj : optimized objective function
  % resL : residuals of the estimation
  % gamma_hat : vector storing the estimates of the intercept and slope
  % coefficients
  
  n=length(y);
  
  %******** MM tolerance and precision variables ****************************

  tn=toler/n;
  e0=-tn/log(tn);
  epsilon=(e0-tn)/(1+log(e0)); %  epsilon is the smoothing parameter.  
  % As suggested in Hunter Lange paper, it
  %  is chosen to roughly satisfy -epsilon(log epsilon)=toler/n using a
  %  Newton-Raphson step above.
  
  %******** Algorithm initialisation *******************
  iteration=0;
  change=realmax;
  res=y-exp(x*gamma); % Compute residuals for current gamma.
  weights=1./(epsilon+abs(res));
  v=1-2*q-res.*weights;
  lastsurr=-sum(res.*v)-sum((1-2*q).*res); % Last value of surrogate (up
                                             % to a constant).
  
  %******** Algorithm loop *******************            
  
  while change > toler 
	iteration = iteration + 1;
    J=x.*(exp(x*gamma)*ones(1,p));  % J is the first differential matrix
       	step = -inv(J'*(weights(:,ones(p,1)).*J)) * (J' * v);
	    gamma=gamma+step;
       	res=y-exp(x*gamma);
        surr = sum(res.^2.*weights) + sum((4*q - 2).*res);
        
        %  Now do the step-halving procedure to ensure a decrease in the
        %  value of the surrogate function.

        while surr > lastsurr
       		step=step/2;
       		gamma=gamma-step;
       		res=y-exp(x*gamma);
       		surr = sum(res.^2.*weights) + sum((4*q - 2).*res);
       	end

       	change=lastsurr-surr;
	    weights = 1./(epsilon+abs(res));
	    v = 1-2*q - res.*weights;
	    lastsurr = -sum(res.*v)-sum((1-2*q).*res);
  end
  gamma_hat=gamma;
  resL=res;
  obj=sum(q.*res)-sum(res(res<0));

end